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Abstract. One way of formulating optical tomography is to use the Born 
series, which is given by the Green’s function for the homogeneous medium. 
The three-dimensional Fn method is a numerical method of obtaining the 
specific intensity of the radiative transport equation expressed as a superposition 
of three-dimensional singular eigenfunctions. In this paper, we compute 
the Green’s function for the radiative transport equation using the three- 
dimensional Fn method without making diffusion approximation. To illustrate 
the present formulation of optical tomography, we consider optical tomography 
with structured illumination. 
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1. Introduction 

We will consider the inverse transport problem that recovers the absorption coefficient 
from boundary data which are directly measured in the spatial-frequency domain. 

For optical tomography reconstructing the absorption coefficient in the half¬ 
space and slab geometries, it is known that plane-wave decomposition is useful 
[ll [22l [54]. Different numerical algorithms to compute the Green’s function for 
the radiative transport equation as a sum of plane waves have been developed 
[ililllllllllSaiMlllIlliaiinililllilllSlISS]. The use of plane waves reduces the 
three-dimensional problem to a one-dimensional problem. For example, in [mini] the 
searchlight problem was considered using Williams’ pseudo-problem approach |54] and 
the Fjsf method [2T1 |48l [50] . The pseudo-problem was directly used for the searchlight 
problem m 156] . By using discrete ordinates and plane-wave decomposition, Kim 
gave the Green’s function as a sum of eigenmodes which are labeled by eigenvalues 
appearing in the corresponding one-dimensional problem [28] . Markel showed that 
eigenmodes in the Green’s function with plane-wave decomposition are efficiently 
computed with the help of spherical-harmonic expansion [41]. Markel’s method is 
called the method of rotated reference frames and has been intensively developed for 
the three-dimensional radiative transport equation in the half-space and slab geometry 
[301 El [31131 [Ml [35]. The recently proposed three-dimensional Fjv method also uses 
plane-wave decomposition |38] . 

In the pseudo-problem by Williams, coefficients of the one-dimensional transport 
equation depend on angular variables. Hence inverse problems with this approach 
are not easy. The discrete-ordinates expression by Kim has been applied to 
optical tomography. In light of three-dimensional singular eigenfunctions [37] . this 
discrete-ordinate method can be viewed as a three-dimensional extension of the 
analytical discrete ordinates developed in one-dimensional linear transport theory 
H El 1 0 0 [H [H. Eigenfunctions of the one-dimensional radiative transport 
equation for isotropic scattering discretized by discrete ordinates are not singular 
due to the interlacing property between eigenvalues and nodes of discrete ordinates 
[53] . However, this interlacing property is lost in three dimensions and the eigenmodes 
become singular depending on the frequency of plane waves. Although the method 
has been numerically proven to work [191 EHl [131 IMl ESI ES] , it ignores the contribution 
from the singular parts of singular eigenfunctions, and the exact value of the specific 
intensity is not obtained even in the limit that the number of nodes for discrete 
ordinates goes to infinity. Optical tomography using the method of rotated reference 
frames was proposed m and is verified by simulation and experiment |40] . Markel’s 
method of rotated reference frames, however, sometimes suffers from numerical 
instability. In particular, the method does not work when the source term has high 
spatial frequencies, which are necessary for reconstruction with high resolution in 
optical tomography with structured illumination [38]. 

Seeking a mathematically satisfactory way to construct transport-based optical 
tomography, we adopt the three-dimensional method |38] in this paper. By using 
the radiative transport equation instead of the diffusion equation, we will generalize 
the noncontact diffuse optical tomography with structured illumination proposed by 
Lukic, Markel, and Schotland [35], which was also experimentally justified [28] . 

In a typical noncontact optical tomography, a source-detector pair consists of 
a point source by raster scanning a collimated laser beam and a pixel of a GGD 
camera. Gompared with the use of large arrays of optical fibers, noncontact optical 
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tomography can readily acquire a large number of source-detector pairs. The Fourier 
transform is performed to the data from such source-detector pairs. By the use of 
spatially modulated beams of structured illumination we can omit the Fourier 
transform for source positions. Thus the number of measurements can be reduced 
compared with raster scanned point sources. In this setup of optical tomography with 
structured illumination, improvement of spatial resolution is shown [7]. Absorbers of 
different structures can be reconstructed by using bi-dimensional source patterns |12) . 
The use of angular-dependence of structured light reflectance [29] and dense sampling 
m were proposed. 

Albeit the use of structured illumination is promising, optical tomography with 
structured illumination has been limited to diffuse optical tomography, which is based 
on the diffusion equation obtained in the asymptotic limit of the radiative transport 
equation. The limitation stems from the fact that calculating the specihc intensity 
of light by solving the radiative transport equation is computationally hard. In this 
paper a reasonably low-cost computation is achieved with an algorithm based on the 
analytical solution to the radiative transport equation. 

Let r = {p,z) be a vector in where p G is a vector in the x-y plane. 
We consider a medium occupying the half-space {z > 0) in which light propagation 
is characterized by the absorption parameter pa and scattering parameter ps- We 
assume that pa depends on r and ps is a positive constant. We write Pa(j) as 


Po(r) = Pa + Spa{r), 

where /io(r), Pa are both positive. We introduce 77(r) as 



where w = ps/pt (0 < -ntj < 1) is the albedo for single scattering. We define 


r 3 = {r e r3. > 0} , = {§ G ±p > 0}, 


where p is the cosine of the polar angle of s, and 


r± = {r G R^, s G S±; z > 0} . 


Let J(r, s) (r G R^, s G S^) be the specific intensity of light at r traveling in direction s. 
We introduce the total attenuation pt = Pa+Ps- In this paper we will take £t = 1/pt to 
be the unit of length. The specihc intensity I obeys the radiative transport equation. 


[ /(r, s) —>• 0, z —>■ oo, 

where /(p, s) is the incident beam. The scattering phase function p(s, s^) is normalized 
as 


1 




z —>■ oo, 


(r, s) G R3 X S^, 



p{s , s) ds =1, s G S^. 


Assuming rotational symmetry we model p(s,s^) as 
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where L > 1, /3o = 1, 0 < /?; < 2Z + 1 for Z > 1, and Pi and Yi^ are Legendre 
polynomials and spherical harmonics, respectively. In the case of the Henyey- 
Greenstein model, we have L —>■ oo, /?; = (21+l)g*, where g is the scattering asymmetry 
parameter. 

Let us consider the Fourier transform of rj: 


^(q, z)= [ z) dp. 


Using fj we introduce 

Vb{p,z) = -77^[ e-^^'PxcienB{<i)v{<l,z)dq, ( 2 ) 

(27r)^ 7 r2 

where is a subdomain in the first Brillouin zone (see ^2.21) and x is tbe characteristic 
function such that y = 1 for q G fls and y = 0 otherwise. We note that rjB- 

As is explained in ^2.31 for each q we have the following linear inverse problem 
within the first Born approximation. 


2^(qo>q) 


K{(\io.z\q)f]B{q,,z)dz, 


Jo 

where qp is the spatial frequency of the incident beam and P is obtained from 
boundary measurements. Thus the quality of reconstruction of q in this inverse 
problem is determined by the kernel K. So far in most research including |35| K 
has been calculated making use of diffusion approximation. In the present paper, we 
will compute K based on the analytical solution to the radiative transport equation, 
without using diffusion approximation. The Green’s function is computed with the 
Fn method, which was recently extended to three dimensions [38] . 


2. Optical tomography with structured illumination 


2.1. Structured illumination 


Let us consider a structured illumination in the half space, i.e., the incoming beam is 
given by 


/(P,s) =/b^(p,s) = Jo[l + AoCOs(qo-p + So)]d(s-So), po e (0,1], 

where Ip is the amplitude, Aq is the modulation depth, and Bq is the phase of the 
source. Since 

^ 2/o(p,s)-(l+i\/3)/_2^/3(p,s)-(l-i\/3)/2V3(p,s) =/q^(p,s), 

where 


/qo(P>s) = /oe - sp), ppG(0,1], 

we will use /(p, s) = (p, s) for the boundary condition [35]. Then we have 

/(q,s)=/' e*'i '’/(p,s)dp = (27r)2/o(5(q-qp)(5(s-sp). (3) 

dR2 

We will reconstruct q(r) by measuring the exitance or hemispheric flux J-|_ defined as 
follows on the boundary. 

J+{P) = / pl{p,0,-s)ds. 

Jsl 
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Let us introduce the data function D{p) as 


D{p) = [ f f pG{p,0,-s;r',s)ds 

JR3 Jsi 


dr'ds, 


(4) 


where I^^\p) is the specific intensity for ? 7 (r) = 0 in and tjb is introduced in ([2]) 
and the Green’s function G(r,s;ro,So) satisfies 

' s • VG(r, s; Tq, Sq) + G(r, s; Tq, Sq) 

= zu p(s, s )G(r,s';ro,So)(is'+ (5(r-ro)(5(s-So), (r,s)eK^xS^, 

G(r,s;ro,So) = 0, (r,s)er_, 

, G(r, s; ro, So)-J> 0, z->• oo. 

Within the first Born approximation (i.e., we assume that |7?| is sufficiently small), the 
data function is given by 

i?(p)« Jf(p)-J+(P). (5) 

We note that the right-hand side of ([5]) is given by the so-called Born series and D{p) in 
O is the first term of the series. Thus D{p) is obtained through the exitance measured 
on the boundary. We can reconstruct jys by solving the linear inverse problem in (j4]). 


2.2. The data function 

The data function (U) is obtained as 


Dip) = Polo [ [ [ mG(p, 0,-s;r',s')ds 

^ G(r',s';p",0,so)e-*‘io'^'"dp' 


dsir') 


dr ds. 


Suppose that there are x Nd detectors on grid points p = {xi, yj ) with spacing hd 
(= Xi+i — Xi = yi+i — yf). We consider the Fourier transform 

D{ci)=Y,D^''’D{p), 


G(p,z,s;p',z',s') = 


(27r)2 




/R2 


G(z,s; z',s';q) dq. 
2 


Noting the Poisson summation formula Xp'^(‘ll + P)) where p 

denotes reciprocal lattice points {2TTi/hd,‘2TTj/hd) (*,j S Z). We assume that TLb 
is in the the first Brillouin zone [—7r/hd, Tr/hJ x [—7r//id, Tr/hJ. We obtain 

Polo 


D{q) = 


hi 


■E 


§2 Jo 


tsl 


^G(0, —s; z', s ; q -1- p) ds 


X ’?B(q + P - (lo, z')G{z', s'; 0, So; qg) dz'ds. 
If q — qQ S Hb, we have 

i)(q) = / [ [ pG*iptz',-s;0,s;q)ds 

^d Jo JS2 [j§2^ 

X G(z', s'; 0, So; qo) ds'dz'. 


( 6 ) 


dB(q- qo,^') 

















We define 
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/W(r,s) = f f G{r,s;r',s)^j,'f^'‘\p',s)S{z')dr'ds\ 

Jsl 

/W(z,s;q)= [ G{z,s;0,s]q)pf^^\q,s)ds, 


with boundary values 

/(^)(q,s) = l, /(^)(q,s) = 5 (s-So). 

We have 

r pOO r ^ 

D{q) = -^ / qQ,z')i'''^\z\s-q^)dsdz'. 

n-d Jo Js‘^ 

We will numerically compute using the three-dimensional Fn method. 

We first subtract the ballistic term: 

/W(r,s) = 4*^(r,s) -b/W(r,s), 
where i = 1,2. Here i]^\y,s) satisfies 


s • (r, s) -b II"’ (r, s) = 0, (r, s) e x 


4*^ (!■>§) = f’'"\r,s), 


(r,s) e r_, 


and Is'^\r,s) satisfies 


■ V/W (r, s) -b /W (r, s)=zu [ p(s, (r, s') ds' + (r, s), 

JS2 

')(r,s) = 0, 


where 


(r, s) =-07 f p{s,s')lj^^\r,s') ds'. 
7s2 


As z —>■ oo we have ^ 0. We note that 

4 * 4 i',s) = e"VW(p,s). 

We expand and with spherical harmonics: 

/’CO CO 






s' G S^, 


m— — co Q:=0 

l=\m\+2c 


CO CO 


s' eEJ_, 


and 


m — — CO OL=0 

l = \Tn\+ 2 a 


/CO CO 






I' g®2 




m— — co Q!=o 

i = |m|+2a 


oo oo 


Y Y 4 m(q 0 )^')^im(-s'): s'gS?_. 


m— — co “=0 

i = |m|+2a 


Note that !;„(—§') = (—!)*¥;„(§'). Since Sq G and — s G in db]), the ballistic 
terms do not contribute to 7i(q). 


(r,s) G X S2, 
(r,s) G r_. 
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2.3. Reconstruction 

By redefining q (q — Qo —q) we have 

poo 

2 ^(qo:q)=/ K{q^a,z;q)fiB{(i,z)dz, 
Jo 


q e 


(7) 


where 


and 


77(qo,q) = i)(q + qo) 


lo' 


oo oo 


Kiqo,z;q)= Y. Y. 

m — — OC O'=0 m' — — oo a' = Q 

l = \m\+2ot l' = \m'\+2ct' 

X km(qO> z)b^l’^,{q + qo, z) + 6|^(qo, + qo. ^)1 [ >7m(s)y;„,(s) ds 


/Si 


oo OO 


= E E E 

m= —oo a:=o q'=o 

l = \m\+2c. l' = \m\+2c.' 


aiJ(qo: (q + qo: + ^iJ(qo> ^)4m (q + qo. 


BJ?,. 


( 8 ) 


Here, 


_ 

Oil/ — 


— m)\(l' — m)\ 

2 V (^ + to)!(/' + to)! Jq 


pr{ti)Pir{fi)diJ. 


Thus the kernel K is obtained using the three-dimensional Fn method. 

Viewing © as a linear matrix-vector equation we can express © as 

\V{q)) = K{q) |i7B(q)), 

where (qo|P(q)) = T>(qo,q), (qol 7t:(q) l^) = K{qQ,z]q), and (2:|i7B(q)) = fiB{q,z). 
We can compute |?7B(q)) with singular value decomposition. We obtain 

VBiq, z) « (z|it:+(q) \V{q)) , 
where is the pseudoinverse such that 

K+ = (kk'i) ^ . 

Here f denotes the Hermitian conjugate and reg means that the pesudoinverse 
is regularized as explained below. Let cr‘j{q) and |nj(q)) be the eigenvalues and 
eigenvectors of the matrix M (q) whose qp-qo element is given by 

poo 

(qol M(q) |qo) = / K{qf),z-q)K*{q'Q,z;q)dz. 

Jo 

If we use the truncated SVD and take only singular values greater than a threshold 
value ctq, ?7 is obtained as 

1 I' -*q-p cr-2(q)^ (r;j(q)|qp)T>(qQ,q) 

Uo 


?7(r) 


<^7 ><^0 


(2^)" Jns 
E-^*(^o.^;q)(qokj(q)) d<i- 


^0 


X 
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Figure [I] shows the reconstruction of a point absorber with 

Tj^r) = r]aVaS{r - Va), (9) 

where it^a = (0,0,2 cm), 14 corresponds to the volume of the absorber, and r]a is a 
positive constant. In the middle panel for z = 2 cm, the absorber placed 2 cm away 
from the surface is reconstructed while the reconstructed ry Ri 0 in the left panel for 
z = 1 cm and the right panel for z = 3 cm. In what follows, we will consider the 
above-outlined reconstruction algorithm with which Fig. [I] is obtained. 




0 . 


2 . » 10 "’® 


Figure 1. Reconstruction of tj for the point absorber placed at {x,y,z) = 
(0,0,2 cm). The panel shows r] on planes parallel to the x-y plane at, from the 
left, z = 1cm, z = 2cm, and z = 3cm. The field of view is 5.1cm X 5.1cm 


3. Preliminaries 

In this section we introduce our notations. We will calculate a|,^(q, z), a|^(q, z), 
and z) to obtain iF in dH) 

In 113.11 we introduce polynomials and p™. In 113.21 and 113.31 Case’s singular- 
eigenfunction approach and singular eigenfunctions in three dimensions are explained, 
respectively. 

3.1. Polynomials 

Definition 3.1. We introduce hi (I = 0,1,...) as 

hi = 21 + 1- Tz^PiX[o.L]{l), 

with X[o,L](0 tf’-G step function (x = 1 for 0 < I < L and x = 0 otherwise). 

Definition 3.2 (Refs. |161 117)1. The normalized Chandrasekhar polynomials 
(m > 0, I > m, V G are given by the three-term recurrence relation 

i^higTiv) = sJ[l + lY-m^gr+M + 
with the initial term 

m, N ^ (2111-1)!! ^ \/(2to)! 

^J(2m)\ 2.^m\ 
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Definition 3.3. The polynomials pY^{p) (m >0, I > m) are introduced as 

pT{t) = (-1) 


' {I — to )! 




{I — to)! d” 


{I + to)! dp' 


:Pi{p), 


{I + to )! 

where Pi{p) is the Legendre polynomial of degree I and Pf'(p) is the associated 
Legendre polynomial of degree I and order to. 

The polynomials satisfy the three-term recurrence relation 

.JP-m^pT_M - {21 + 1)ppT{p) + ^{l + iy-m^pT+M = 0 - 

3.2. Singular eigenfunctions for one dimension 

In one-dimensional transport theory, singular eigenfunctions (j)"'{i',p) are given by 


^ N = —P ^ ^ ^ ('^)(1-m) S{iy-p), 

where V denotes the Cauchy principal value and 

L 

g^{v,p)= Y, Pipr{p)9ri^)- 

l—\m\ 

Here |to| < L and € K are eigenvalues; z/ has discrete values > 1, 

j = 0,1,— 1) and the continuous spectrum between —1 and 1. The number 
of discrete eigenvalues depends on w and j3i. The function A™(z/) is given by 


2 7_i V - p 

Singular eigenfunctions are normalized as 


{l-p‘^)\'^\dp. 


y'%-(z.,Az)(l-Az")'™' dp=l. 


We note that 


9r{^) = (- 




Discrete eigenvalues are roots of A™, i.e., A'"(z/™) = 0 (w G C), where 

2 w-p 

We have the following orthogonality relations 0112145] 

J^^p(r{T^,T)(r{^',p)dp=M'^{v)5,u', 

where the Kronecker delta 5^^' is replaced by the Dirac delta 5{v — v') if z^, v' are in 
the continuous spectrum. The normalization factor Af'"'{v) is given by 


V'"(z/) = 


dw 


zzA™+(zz)A'"-(z.)(1-z.2)-I'"I, 

where K"'^{v) = lim£^Q+ K"'{v ± ie). 

Finally we introduce 


V G (-1,1), 


<^^{s) = <r{v,T){l-9^) 


| m |/2 
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3.3. Singular eigenfunctions for three dimensions 

Let '(/'(§) G C be a function of s G By the operator TZ^ defined in [35], angles in 
TZ^ ipis) are measured in the rotated reference frame whose z-axis lies in the direction 
of a unit vector k G (k • k = 1). 

If '(/’(§) G C has the form 

OO I 

m^lm (s)) i’l m G C, 

l—O m— — l 

then we have [11123111] 

OO I I 

^k^(§) = E E E 

l—O rn— — l m' ——l 

OO I I 

^£'V'(§)= E E E 

m— — l m' ——l 


where 9^ and are the polar and azimuthal angles of k in the laboratory frame, and 
Wigner’s d-matrices. 

We note that TZ^ §•§' = §•§' and TZ^ ^ = s • k. 

Definition 3.4 (Plane wave decomposition). Complex unit vectors k(z/,q) G 
G M, q G are given by 


k{v, q) 


Jvq \ 
kz{vq) ) ’ 


where 9 = |q| and 

kz{yq) = Vl + {vqY- 

It turns out that q)) functions of vq. We write 

dLm'[iT{vq)\ = C^Lm'(%(;.,q))- 

These d-matrices are computed using recurrence relations [Ml [33138]. To calculate 
square roots such that 0 < arg(i/i’) < tt for all z G C [Ml [33] , 

We have 


‘/’k(!^,q) 


yjq, for V > 0, 

(fq + TT, for V < 0, 


where (pq is the polar angle of q. 

Three-dimensional singular eigenfunctions are obtained as 


^ i/ — s • k 




s-k). 


where k = k(i/, q). 

The following orthogonality relation holds. 






'(§ 




*(§)) ds = 2'nkz{vq)M{v)5vv'5r, 
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Since the general solution is given by the sum of a particular solution and a linear 
combination of eigenmodes, we can write the Green’s function as [5] 

G(r, s; ro, § 0 ) = Gfree(r, §; Tq, Sq) 


c/IR2 ^ — _t 1 o—n 


L (M^-l 


m— — L K. j—0 


+ I q) dy \ dq, 


where Gfree(r, §; ro, § 0 ) is the free-space Green’s function and A”^{v) are 

some constants. We can readily check that the above expression satisfies the radiative 
transport equation and the boundary conditions. We note that the free-space Green’s 
function is obtained as m 

Gfree(r,s;ro,So) = 7 ;^ / e“*‘i'(^^^«)Gfree(z,s;zo,So;q)rfq, 

(27r)^ Jr2 

where 


L (M^-l 


1 






m— — L K j=0 


+ 


■7^c 




Jo 2Trkz{vq)JJ'{i^) 

Upper signs are chosen for z > zq and lower signs are chosen for z < zq. 


4. The three-dimensional Fpf method in the half space 


Let us consider a homogeneous medium (? 7 (r) = 0). Let s) denote the specific 

intensity in the medium. We can split 7^°^ into the ballistic term If, and the scattering 
term Is'. 


=7f, + 7,. 


Here, 


and 


We obtain 


( § • V/b -|- 7b — 0, 

< 4(r,s) = /(p,s), 
[ 7b(r, s) 0 


(r,s) e R3 X S2, 
(r,s) € r_, 
as z —>■ 00 , 


s • V7,,-7 7,, = rn f p(s, s')7,,(r, s') ds'-7 S'[/](r, s), 

Is{r,s) = 0, 

7„(r, s) -J> 0 


(r, s) G X S^, 

(r,s) G r_, 
as z —>■ 00 . 


7b(r,s)= / S{p-p')- - - — d{s-s)pf{p',s)d{z')dr'ds=e ^/(p, 

dRixSi M 
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for s G , and /{, = 0 for s G . Hence 


5'[/](r, s) = tz7 / p{s,s')Ib{r,s')ds' =we ^ [ p{s,s')f{p,s')ds'. 
J§2 Jsl 

The scattering term is given by 

Is{r,s) = f G{r,s;r',s')S[f]{r',s')dr'ds'. 

x§2 


We consider the Fourier transform h of h- 


isiq,z,s)= / e*'i'^/s(r,s)dp 

- ^ ^ G{z,s]z',s\a^S[j]{(\,z',s)dsdz' 


'0 JS^ 


G{ree{z, s; z', s'; q)S[f] (q, z', s') dsdz' 


10 dS2 




m— — L \ i—0 


with some coefficients A^^iy). 

Let us introduce the notation 

V L(''V-|m|)/2j 

E= E E ^ 


Im m— — N 0 

Z=|m|+2Q 


where N {> L) is the maximum degree of spherical harmonics in the expansions. For 
s G S^, we expand Ig as 

is(q,0,-s) « ^ cim{q)Yim{s), 

Im 

Is{q-,z,-s) K. ^ 6;„(q,z)ri„(s), 

Im 

Is{q,z,s) « ^ a;m(q, z)rim(s). 

Im 

For simplicity we set 

N = L. 


Definition 4.1. We define (—L <m<L) as 


xm _ 

?.7 “ 


TT 

cos I — 


j - M”" + 1 
2 TV™ - M™ + 1 

row ' 


where 


NZ 


N — \m\ 
2 


+ 1 . 


We drop the superscript m if there is no confusion. 
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Let m' be an integer between —L and L. We introduce 


/S2 


= X[-iA + (GO- 


We have 



o 

s) ds 

= if [/](-?., q,o), 


( 10 ) 


o' 

)K-^ 

s) ds 



( 11 ) 

J §2 V y 

) Isi<l,Z, 

s) ds 

= if[/](-?.,q,^), 


( 12 ) 

/ A^f^k <*(§)' 

J §2 V / 

) 4(q,^, 

s) ds 



-kzi^jq)z/ij 




+ if [/](^j>q:^), 


( 13 ) 


where 

Lf[f]{-i„ci,z) = 


and 

Lt[mj,ci,z) = 


n 


W^j 


■k(-C,.q) 

L I 




"^[/Kq, z',s )dsdz' 








/o Js^ 




— (^-kz{^jq)z/^j _ p- 


1 


w^j 


L I 


H H / /(q>s)Li 






In particular for 


/(q,s) = <^ ’ 

( 5 (s-z), 


we have 


LT'[f]{-^j,<i,z) = e-^ ^ (-l)'AdL'[*T(6<?)]ffr'te)x 

2fe- + 


and 


Lf[mj,^,z) = - e-^) 

27r f Pi{fi)dn, 

Jo 


W^j 


2(Ci kzi^jq)) 


l^idom'[>-Ti^j<l)]9r' i^j) 


3 ^'' l—\m'\ 


*m{s)ds. 
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where we used Yi^iz) = {21 + l)/(47r)i5mo- Here, 

1 , 

0 , 

^0 d! 


[ Pi{fi)dn = 

Jo 


i-lf 


-l)/2_ 


1{1 + 1){1-1)\V 


I = 0 , 

even I (^0), 
odd 1. 


By noticing = (—we see that L™ [/] and L™ [/] are independent 

of ipq and independent of the sign of m'. 

From (fTOl) we obtain 


-(-ir' / (§))/,(q,0,-s)ds = Lr'[/](-G,q,0), 

J 

where we used (s) = (—1)™ (“§)■ We introduce 

Explicit expressions of are found as follows [38]. 

= k^{^3d)sj^^{-^)'^d^rnn.AiT{£.jq)\ (^\/(^ + 1 )^ “ ( 0 ) 

’ m" — — l 

dm",m'-i {I - m"){l - m')5z-\(G) “ + m' + 1)(? + wOffEfe)) 


Im 


-I- 


+ 6m",m' + l (^V{l-m' + 1){1 - TO')5iTl(6) - \/G + to")(^ + w')5i-l(0)) 
_ yl-)!™' 

/m ’ 

where we used g’^{—v) = (—and 




47r (1 + to)! 


( 2 | to '| - 1 ) 1 ! 




? 7 i" = —|m'| 

5™ costp) 

+ kz{ijq)^^ - cost/? 


/(I™' 

— to")! 

J (\m' 

1 +to")! 




When numerically evaluating the above integral over s, we use the Gauss-Legendre 
quadrature for g and trapezoid rule for ip. We note that 

/_-I \m „(±)jm' 

'^L-m ~ V '^Im 

For each , q, we have 




Im 
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We have 




Im 


Im 

Im 

Im 

The above relation implies 

Qm(q) = Cl-miq) = (“ l)™Qm (<?) • 

That is, 

N 

E E + (1 - SmoK-irjlt^r] cimig) = -Lr'[/](-c„q,o). 

iTi—0/^|m|,|m|+2,... 

We obtain from m and (US 

fi (7^£ <i>?^'*(§)) /s(q, §) ds = ^ (Te^ $^'*(§)) /.(q, 0, s) ds 

+ Lfm„q,z). (14) 

Together with (fT^ . we have 


Im 


Im 




Im 


Im 


= +L^'[/](e„q,^), 


where we used Yim{—s) = (—l)*Yim(s). Let us define 

vr' = -L-'[/](-e,.,q,z), 






and 


We can show that and are independent of (p^ and satisfy v(’ ™ = v("^ and 
/ / 

^2 = ^ 2 ™ ■ Similarly to Cim(q) we have 

atm(q, z) = aim{q, ai-m{q, z) = {-l)"'aim{q, z), 

blm{q, z) = blmiq, -2)6“*’”'^'’, k^^jniq, z) = (-l)™ 6 im (g, z). 

= (-I)'J,+ (1 - 

+ (1 - s.M-irjYY ■ 


We write 
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We obtain 


... 

3 

1 

. . . ... 

••• ••• 

jTJi',Lm 

••• ••• 

jm',1771 


\ / 


aimiq.z) 


bim{q,z) 


V 


( : \ 


/ 






V : / 


where 0 < to' < L, 0 < j < IV^w — 1, 0 <to<-/V, Z = to + 2a, 0 < a < | |. 

In [Appendix A[ the three-dimensional fjv method for the slab geometry is 
explained. 

Finally we can rewrite (H)) as 

N l{N-m)/2\ L(JV-m)/2j 

771 — 0 “ = 0 q' = 0 

i = |m|-f2Q' l' = \m.\-^2a' 

X z)b''P ^(|q + Qol, z) + 6|^(<7o, 2)4m*(1^ + ^oU) 

X [^mO + 2(1 — Smo) cos (</5q+qg — V^q^)] 


5. Simulation 


Let us test our optical tomography developed in lJ5]using the following simple example. 
In this section we describe how Fig.[T]is obtained. We suppose diffusion approximation 
works well. Although the diffusion equation is obtained as an approximation to 
the radiative transport equation, we will calculate 2 ?(qQ,q) by solving the forward 
problem of the diffusion equation to avoid inverse crime. [Appendix B| is devoted to 
the computation of the data function 2 ?(qQ,q). 

We consider a random medium with A = 3, 

/io = 0.05 cm“^, = I00cm“^, 5 = 0.9005, (15) 


for which we have £* = 1 mm. Here £* is the transport mean free path given by 
£* = {^J-t - ^^s9Y ■ 

We assume that a point absorber is embedded at Xa = ya = 0, Za = 20 mm. We set 
rya = 0.0015, Va = 10“®mm3. 


We note that 


4 /v Ma(r)dr //ia ~ 4. We place detectors on a square lattice of 


spacing hd- Detector positions are specified by 


We set 


ixi,yj) = iihd,jhd), i,j = -Nd,...,Nd. 


hd = l mm, Nd = 25, 

for which the field of view becomes 51mm x 51mm. Let q(y)(i) denote the x- 

and ^-components of q. We put W = [2{Nd + Ng) + l]hd: and 


q 




2tt 




hj = -Nd, ■■■,Nd. 
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For Qq we use 


We set 


„( a :)(0 

% 


2tt 


lv)U) 


27T . 


% VF'^’ 


hj 


= -N. 




Ns = 10 . 

Note that ~ O.lmm”^. In this demonstration we chose N = L = 3. 

The first 10 largest singular values are typically used for the truncated SVD. Finally 
we put Iq = 1, §0 = z. 


6. Concluding remarks 

Optical tomography based on the radiative transport equation was considered using 
the three-dimensional FN method as a precise and stable numerical algorithm. It 
is straightforward to extend the Fjv method for the half-space to the slab as is 
described in [Appendix A[ Thus optical tomography in the slab geometry can be 
similarly formulated. 

Although the first Born approximation is employed in this paper (see ([5])), recently 
optical tomography with higher-order nonlinear terms has been developed by the use 
of a recursion algorithm for the inverse Born series of the radiative transport equation 
[3^ . The present scheme for calculating the kernel K with the three-dimensional 
Fjv method can readily be extended to such nonlinear inverse problems, which is an 
interesting future problem. 


Appendix A. The Fjv method in the slab geometry 

We consider a slab of width The radiative transport equation is given by 

s • V/(r, s)-|-/(r, s) = 117 / p{s,s)I{r,s)ds, (r, s) G x 

/(r,s) = /i(p,s), z = 0, p G § G §2 ^ 

/(r,s) = / 2 (p, s), Z = Zmax, P G R^, § G §?., 

with some boundary values /i(p, s), / 2 (p, s). The specific intensity is given by 

1 


(A.l) 


'M^-1 


/(r,s) = 


(27r)2 




/R2 


■'E E 




$"*(§)£ 




m— — L \ j—0 


M"”-l 

i-o 

with some coefficients A"*(±z/™), Using the orthogonality relations we have 

'' p (7^fc $-*(§)) /(p, 0, s) ds = (A.2) 


/S2 


P (7^fc $-|(s)) /(p, z^ax, §) ds = (A.3) 

p {TZ^ $"*(§)) /(p, 0, s) ds = 2TTkkq)Afk)A^{0, (A.4) 


diy 
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/ J(p,z^ax,§)rf§ = 

J§2 

[ fx (7^fc$™|(s)) Iip,z,s)ds = 27ri(Cg)Ar(-C)A™(-C)e^^««)^/«, 
Js^ 

[ M (7^k ‘J>r*(§)) §) . 

JS2 

Let us consider the Fourier transform of the specific intensity: 
7(q,z,s)= [ PI{r,s) dp. 


/R2 


We express /(q, s) (s G S^) as 
/(q,0,-s) 


N L(^-|m|)/2j 

^ ^ ar(q)r,™(s), 


n=-N “=0 

Z = jmj+2Q 


7(q, ^max,s) 


/(q,z,-s) 


/(q, 2 :,s) 


N L(JV-|m|)/2j 

E E ^r(q)i 7 n^(§), 


n=-N “=0 

i = |m|+2Q 


N L(^-|m|)/2j 

E E c["(q,z)ri„(s), 


m— — N «=o 

Z = jmj+2a 


W L(JV-|m|)/2j 

E E dTi<l:Z)YUs), 


n^-N «=o 

Z = jmj+2a 


where s G S^. From (IA.2I) and (IA.3I) we obtain 

-i-ir [ /r(7^£<i>4™*(s))7(q,0,-s)ds+ / /r(7^£ $“!(§)) A(q,s)ds 

p (7^£ -^(q> -ZmaxjS) ds 


— g (^Q')2^max/^ 


(A.5) 

(A.6) 

(A.7) 


_(_l)mg-L(««)^max/« f ^(7^- $™*(s)) / 2 (q,s)ds, 

Jsl 

where we used = (—1)'"$™(—s). Similarly from (jA.4ll and (|A.5p we obtain 

f ^ (7^j^$^*(s)) 7(q,0max,s)ds - (-1)™ / /r (7^i^ $“!(§))/2(q,s)ds 

= e-fcA«9)^max/« f ^(7^-$™*(s))/l(q,s)ds 
dS2 

_(_l)mg-fc.(^9)z„ax/« f ^(7^- $™*(s))7(q,0,-s)ds. 

Jsl 

They are rewritten as 

E^i™'(e)ar' + = ET2, 

Im' Im' 


^21^ 


(A.8) 

(A.9) 


Im' 


Im' 
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where 


AZ'i-0= [ (Tet (§)) (s) ds, 

Jsl 

Bz,{o = i-ir f /r(7^ii<i>^(§))^™'(§)d§, 

Jsl 

ET2 = [ M(%<5-5(§))/i(q,§)rf§ + (-ire-"^^«^^"““/« / fi{ni^^f*{s))f2{q,s)ds, 
Jsl Jsl 

E^, = f M(7et$!!*5*(§))/2(q,§)ds + (-ire-^^(««)^““/« / M(7^k$^(§))/l(q-§)'^§• 

Jsl 

We obtain by using (jA.2ll and (|A.6p 

[ M (®)) ^(q’0,s)ds = ^ (%$™j(s)) J(q,z,s)ds, 

JS2 JS2 

from (lA.Sp and (lAjp 

[ = [ /r(7^£$^*(s))/(q,z,s)ds, 

J§2 JS2 

by using (lA.Sp and (lA.Bp 

[ = [ /i (7^£ $“!(§))/(q,z„,ax,§)d§, (A. 10) 

JS2 

and from (jA.4ll and (jA.7l) 

f A^(7^ii$^*(s))/(q,^,s)ds = e-^^(^^)"/^ / /.(7^fc$^*(s)) J(q,0,s)ds. (A.ll) 

JS2 JS2 

We obtain a™, from the linear system (IA.8I) . (IA.9I) . Then c[", d™ are computed 
using (lA.lOl) . (lA.llI) . We note that (lA.lOl) corresponds to (fT^ and (lA.llI) corresponds 
to (fT4|) . 

Appendix B. Forward problem 

We will calculate the data function using diffusion approximation. The unit of length 
is it ■ With diffusion approximation we can express the Green’s function as |43] 

G(r,s;r',s') ~ ^^^(1 + Ts • V)(l - Ts'• (r', r), 

Att 

where = w{l — g) and i* = 1/g* = 1/(1 — wg). The Green’s function is further 
approximated by 

G{z, s; z\ s'; q) ~ ■ q)(l + iGs' • q)G(°®)( 2 ;, z'; q). 


47r 


Here 


and 


G(°E)(r,r') = ^ / 

(27r)^ Jr2 

I -DoV2G(°®)(r, r') + aoG(°^)(r, r') = d(r - r'), z > 0, 
I G - • VG = 0, z = 0. 
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where 


D 


^0 — C/Xq,, 


with c the speed of light in the medium. We choose the extrapolation distance 
{0 < £ < oo) as 


2 


We have |42| 


where 




2DoQiq) 


„-Q(g)|z-z'| _ -Q(q)\z+z'\ 

1 + Q{q)e 


Q{q) = \Jkl + q^, ko = J^ = = \/3(l - w){l - wg). 


We have ko ~ 0.12 mm ^ for optical parameters in (I15|) . In particular, 


G(°^)(z,0;g) = 


p-Q(9)2 


z>0. 


Dq{ 1 + Q{q)£) 

Using the radiative transport equation, the energy density (r) is obtained as 

u(o)(r) = i / /(o)(r,s)ds. 
c 7 s2 

Its Fourier transform has the form 




1 


S2 Jo 


G{z, s; z', s'; q)At'/(q, s')(5(2:') dz'ds'ds, 


where we used Sq = z and ([3|). By diffusion approximation we have 

{2(°^(q,z) ~ (27r)^/o/i'r(5(q-qo)G(°’^^(2;,0;qo). 

The Fourier transform of the following is given by the right-hand side of the above 
equation. 

u(o)(r) = f r')S'(r') dr', 

Jr^ 

with the source term 

S{r) = Ioki'ste-^^°-PS{z). 

In diffusion approximation M(r) = (1/c) Jga J(r, s)(is, where J(r,s) is the solution to 
o, satisfies 

f -iJoV^u(r) -h (ao + cf]{r))u{r) = S{r), z > 0, 

1 u — iz ■ Vu = 0, z = 0. 


We have 


i(r) = f G^°^^(r,r') [S'(r') — cr7(r')u(r')] dr' 

Jr^ 

= tt^'’^(r) —c f G^^^^ (r, r')?7(r')M(r') dr' 
Jr3 

= ■u(‘’^(r) - 7oG(°®)(r,ra)u(ra), 
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where 

7o = crjaVa- 
The Born series is given by 

■u(r) = - 7oG^°®Hr,ra)u^°Hra) 

+ 7o(r, (r,, ra)w(°) (r,) 

- 73g(°’^)( r, ra)G(°^)(r^, ra)G(°^) {va, ra)u^°'>{va) + ■ ■ ■ 

= u(°)(r) -7rG(°®Hr,ra)u(°)(ra), 

where we write m 

__To_ 

1+ 7oG(°®)(ra,ra)' 

Introducing the ultraviolet cutoff (fcgAc <C 1), we define the Green’s function 
G(DE)(r^^ro) as [H] 


G(°^)(r„,r,) = 


1 


47rGo 

1 


C^tt/Ac 


■ dq + 


q Q(gK-l ^-20(,).„ 


0 Q(,<i) Jo Q{q) Q{<i)^ +1 


27r\ 


^ — 2koZa 2 


47rGo 

where Ei is the exponential integral defined by 

roo ^ — t 

Ei(z) = 


Q\^]-ko + - 2Zaiko + y) 


Ac; 


2^0 


coo g-t 

,0) = / -dt 


In numerical calculation we set 
Ac = Vj/^. 

Within diffusion approximation the hemispheric flux is written as 


jr\p)=^ 


r27z pO 




I J ^ [cu(p, 0) — 3Z?oS • Vu(p, 0)] 

2\2 ci 

We then have 

^PE)(0)(^) _ ^|DE)(^) ^ C ^ ^GDE(p,0,ra)w(°^(ra). 

The data function is calculated as 

lo 




7Pb)(0)(^)_^Pe)(^) 


p 

Qh^'Jrllsi -iq„-p„ ^z(q+q„)-p 6 

2(27r)2cr 
ISnqaVatJ.'si'^ 


E' 


-Q(9o)Za 


p-*q' (p-Pa) . 


.gi^-Pa 


1 + Q{qo)( Jv? 

g-Q(9o)^a g-Qdq+qoD^^a 


-Qiq')za 

1 + Wy 


47r£* + SrjaVaCo 1 + Q{qo)£ 1 + Q(|q + QqI)^’ 


where 


Co = <3 ( ^ 1 - ^0 + 


For large Za, we have 

Co ^ <3 


-2kQZa 2 


22a 


- 2zaiko + -) 


- kn. 


dq' 
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